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ABSTRACT 


This revised user's manual describes the details of a general-purpose 
computer program VISCEL (VISCoELastic analysis) which has been devel- 
oped for the analysis of equilibrium problems of linear thermoviscoelastic 
structures. The program, an extension of the linear equilibrium problem 
solver ELAS, is an updated and extended version of its earlier form (written 
in FORTRAN II for the IBM 7094 computer). A synchronized material prop- 
erty concept utilizing incremental time steps and the finite element matrix 
displacement approach has been adopted for the current analysis. Resulting 
recursive equations incorporating memory of material properties are solved 
at the end of each time step of the general step-by-step procedure in the 
time domain. A special option enables employment of constant time steps in 
the logarithmic scale, thereby reducing computational efforts resulting from 
accumulative material memory effects. A wide variety of structures with 
elastic or viscoelastic material properties can be analyzed by VISCEL. 

The program is written in FORTRAN V language for the UNI VAC 1108 
computer operating under the EXEC 8 system. Dynamic storage allocation 
is automatically effected by the program, and the user may request up to 
195K core memory in a 260K UNIVAC 1108/EXEC 8 machine. The physical 
program VISCEL, consisting of about 7200 instructions, has four distinct 
links (segments), and the compiled program occupies a maximum of about 
11700 words decimal of core storage. VISCEL is stored on magnetic tape, 
and is available from the Computer Software Management and Information 
Center (COSMIC). 
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VISCEL- A GENERAL-PURPOSE COMPUTER PROGRAM FOR 


ANALYSIS OF LINEAR VISCOELASTIC STRUCTURES 

USER'S MANUAL 


I. INTRODUCTION 

The general-purpose digital computer program VISCEL is capable of 
solving equilibrium problems associated with one-, two-, or three- 
dimensional linear viscoelastic structures (Fig. A-l). Since the program is 
an extension of the linear equilibrium problem solver ELAS (Ref. 1), its sol- 
ution at the beginning of the initial time step yields elastic solution of struc- 
tures. Basic inputs of VISCEL, thus, are the same as in ELAS; additional 
inputs are, however, necessary for VISCEL, which represent changes in 
material properties and loading in the time domain. Other important features 
of the program include dynamic memory allocation, optional node relabelling 
scheme, boundary condition imposition during assembly of the stiffness 
matrix and its storage within a variable bandwidth. The program is further 
divided into four distinct links, namely, input, generation, deflection, and 
stress links. 

This user's manual describes the numerical problem formulation, input 
preparation, output description, and other relevant details of the program. 

The physical program is available from COSMIC. * 

Volume II of this report is the program manual, which contains the 
lists of variables, subroutines, and flow charts as well as other pertinent 
program information (Ref. 2). 


Computer Software Management and Information Center, Computer Center, 
University of Georgia, Athens, Georgia 30601, telephone: (404) 452-3265. 
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II. 


BASIC CAPABILITIES OF THE VISCEL PROGRAM 


The basic capabilities and the initial inputs of VISCEL are the same as 
the linear equilibrium problem solver ELAS (Ref. 1). In order to achieve a 
self-contained report, this report includes several tables and figures from 
Ref. 1; they are provided in the Appendix in appropriate order. Thus, 

Tables A-l and A-2 describe, respectively, the various structures that can 
be solved by VISCEL and their compatible combinations. Also, information 
regarding various available finite elements is given in Tables A- 3, A-4, and 
A-5. Further, the usual conventions for ordering of element nodes are 
explained in Table A-6. VISCEL can handle any material, namely, isotropic, 
orthotropic, or anisotropic; their input requirements are described in 
Fig. A-2. 

IH. NUMERICAL FORMULATION OF THE LINEAR THERMOVISCO- 
ELASTIC PROBLEM 

Reference 3 gives complete derivation of the numerical formulations of 
the linear thermoviscoelastic problem, whereas Ref. 4 presents details of 
the finite element technique. However, such formulations are summarized 
below in a simplified manner for completeness of this report. 

A. Basic Approach 

The fundamental equilibrium problem in structural analysis can be 
formulated as differential equations with appropriate boundary conditions; 
alternatively, an equivalent extremum formulation may be developed based 
on the principle of minimum potential energy and its complement (Ref. 5). 

In this work, structural discretization is achieved by the finite element dis- 
placement matrix technique, a variant of the well-known Ritz method for the 
minimization of the total potential energy functional i)j associated with admis- 
sible displacement trial functions. The admissible functions are restricted 
to be sufficiently smooth, usually being algebraic or trigonometric polyno- 
mials, and, furthermore, they are required to satisfy essential boundary 
conditions arising from the requirement of geometric compatibility. This is 
achieved by expressing the trial solution in terms of a set of linearly inde- 
pendent known functions and undetermined parameters and then minimizing 
the functional with respect to such parameters. 
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In the finite element method, a structure is discretized by any suitable 
random mesh, and a family of piecewise continuous displacement fields is 
prescribed for each element, which are finally expressed in terms of their 
nodal function values. Such nodal displacements are the undetermined 
parameters to be determined from the extremum principle, the fundamental 
assumption in the procedure being that the total potential energy of the entire 
structure is equal to the sum of potential energies of the individual elements 
(Ref. 4). Such an assumption is valid provided the displacement functions 
and their derivatives of order one less than the highest one appearing in the 
functional are continuous at interelement boundaries; this ensures that values 
of highest derivatives occurring in the total potential energy functional ip 
remains finite (Ref. 4). Obviously, the greater the number of chosen unde- 
termined parameters, i. e. , finer the finite element mesh, the lower the 
value of the total potential energy would be, yielding even better approxima- 
tions. Whereas ip approaches its minimum value from above, the corres- 
ponding strain energy value is always underestimated, and hence the present 
approach computes lower bounds of associated displacements. The finite 
element procedure thus gives a stationary value of ip for the variations of the 
unknown nodal displacements. Because of its resemblance to the piecewise 
Ritz procedure, any particular nodal parameter is only influenced by its 
adjacent elements, and hence the final stiffness matrix is highly banded in 
nature for most practical problems. It can further be shown that the minim- 
ization process of total potential energy of the entire structure with respect 
to each unknown nodal displacement is equivalent to the appropriate summa- 
tion of such process for all individual elements with respect to their nodal 
parameters. For quadratic functionals, the piecewise Ritz procedure for 
each element yields symmetric linear equations in the element displacement 
vector. The minimization process for the entire structure then leads to the 
set of linear, simultaneous equations: 

^K q< r^fq e ^p'=» (1) 
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where 


k® = element stiffness matrix 
6 

q = element nodal displacement vector 
0 

p = equivalent nodal load vector 

with appropriate summation over all elements based on nodal connectivity. 
Such equations are further positive definite for stable structures and may be 
solved by standard processes to yield the undetermined nodal displacements. 
Computation of stresses, etc. , is performed next by the usual procedure. 
The program starts with the computation of element stiffness matrices 
already derived above. 

In viscoelasticity, the creep strain rate or the relaxation stress 
response is dependent not only on the current stress and strain state, but 
also on the entire history of its development in the time domain. Associated 
numerical computation procedures usually adopt a step-by-step incremental 
process, which normally requires knowledge of stress and strain at all pre- 
ceding intervals. This enables computation of stresses/strains at a given 
time implied by some relevant law of the characteristic functions. Usually 
such material properties are strongly dependent on time and temperature. 
The viscoelastic equations are developed as finite difference equations in 
time and finite element matrix equations in space (Ref. 3). This computer 
program is based on linear thermoviscoelastic formulations utilizing a 
11 synchronized" material property concept for thermorheologically simple 
materials. The fundamental assumptions may be summarized as follows: 

(1) Material properties . The material properties may be 

temperature-dependent and are assumed to behave in a thermo- 
rheologically simple way; thus, for temperature changes, the 
characteristic functions, both creep and relaxation, show pure 
shift when they are plotted against the logarithm of time. 

Such materials are better suited for a complete characterization 
over a large range of time and temperature scale since their 
rheological behavior can be described for the entire temperature 
range as a single function of reduced time and temperature. 
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Thus, when any characteristic function, such as the relaxation 
modulus, is plotted against reduced time, all curves will fall on 
the single curve for initial temperature Tg. Hence it is then 
necessary to determine relaxation/creep functions for one tem- 
perature only. 

The shift functions may sometimes be dependent on stresses, 
requiring determination of the shift function at the end of each 
time step. However, such considerations are excluded in the 
present version of the program. The material can be isotropic, 
orthotropic, or general (Fig. A-2), provided they are properly 
defined by experimental results. For this analysis, it is 
required to have a knowledge of the modulus functions (relaxation- 
type functions). Furthermore, the material is assumed to be at 
least slightly compressible. 

The concept of synchronized material properties is that all 
material properties are functions of only one parameter i. The 
same concept applies to external loadings, both mechanical 
and/or thermal. The parameter £ may be time, reduced time, 
or any other suitable variable. Material and load data are con- 
sidered in functional form (Fig. 1), which are to be presented at 
each time step, in the shape of predetermined tabulated values 
obtained either experimentally or derived from analytical con- 
siderations; any interaction between them is assumed to be 
included in such values. 

(2) Linear viscoelastic behavior . Strains are linear functions of 
stresses, but are strongly dependent on loading history, implying 
that if all loads are doubled, all deformations will be doubled too. 
Thus, creep/relaxation laws are linear in stress /strain and as 
such the principle of superposition is valid for such cases. 

Further geometric nonlinearies, e. g. , large strain or large 
deformations, are not considered for the current analysis. 

(3) Deflection boundary conditions . Deflection boundary conditions 
are assumed to remain unaltered throughout the entire time 
domain of computation, being fixed initially at the beginning of 
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the initial time step. The solution at the beginning of such initial 
time step corresponds to the usual linear elastic analysis of the 
structure. 

B. Analysis Review 

Numerical formulation of the step-by-step linear thermoviscoelastic 
analysis procedure for quasi- static problems may now be summarized. A 
basic assumption in the analysis is that the materials are thermorheologi- 
cally simple in nature. Such an assumption is necessary so that the charac- 
teristic functions may be singly defined for the entire temperature range in 
the time domain. The usual field equations for viscoelastic materials may 
then be extended for the thermoviscoelastic case. This is achieved by intro- 
ducing the concept of a "reduced time" when all characteristic functions ful- 
fill the same time-temperature shift and can be represented as a function of 
reduced time: 


/*t 

« (x h"l = a[T(x“7j] (2 > 

in which a(T) is the time shift function usually determined experimentally as 
a function of temperature T only. Such shift function dependence on time t 
and position x^ within the material region is implicit through T, and may be 
sometimes described by the well-known Williams-Landel-Ferry (WLF) 
equation. Relationship (2) signifies that all the characteristic functions, 
such as relaxation moduli of a thermoviscoelastic material at any arbitrary 
temperature T corresponding to time t, may now be expressed by their 
behavior at reference temperature Tq on the new reduced time scale £. 

Each relaxation modulus, signifying relaxation stress variation for unit 
strain applied initially, may then be expressed as 

E 5h<‘> = E ij^> (i = j = k = i = 1, 2. 3) (3) 

Eijk^(t) being the general anisotropic relaxation moduli having 21 independent 
components. The constitutive equations for the usual viscoelastic case may 
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be derived by approximating strain variations by the sum of a series of 
step functions, which corresponds to a series of relaxation displacement 
inputs. Constitutive equations are obtained, from superposition principles, 
in the form of hereditary integrals. For the present thermoviscoelastic 
case, the constitutive equations may simply be derived from such relations 
for the corresponding viscoelastic case by utilizing Eq. (3): 


o\.(x, ,t) = 
ij h’ 




X 




e kf (x h’ T) “ *k (x h’ T)0(x h' T) _ 


(4) 


which may be rewritten as 


°"ij E ijk/^ e k£(0) + / E ijkf ( ^ " dr (e k^ ' a k£ 0) dT (5) 

J 0 

in which is the initially induced step strain at t = 0, the correspond- 

ing first tefm being the effect of such initial strain at time § (x^, t). The 
kernel of the hereditary integral E.^fKx^ t) - |(x h , T )] may be considered 
as the memory function transforming the influence of pulse strain at time r 
to the time instant t. In addition to Eq. (4), two more equations are 
required to completely define the field equations; 

(1) Equilibrium equations 


<r. . . + f. = 0 
1 J > J i 


(2) Strain -displacement equations 


( 6 ) 


e. . 



+ u. .) 
h i 


(7) 


where T is the body force component per unit volume. The field equations 
may next be expressed as incremental field equations when the time domain 
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is subdivided into arbitrary intervals At(m). Equations (6) and (7) and the 
stresses of Eq. (5) then take the following form: 


A<r. . . + Af. . = 0 

ij(m) , j i(m) 


Ae../ . = rr (Au., . . + Au., , .) 

ij(m) 2 i(m),j 


( 8 ) 


°"ij ( n ) 


" E ijk^ ( ^(n) ) e kf(0) 


(n) 


E ijkf ( ^(n) 


- (9) 


Finally, Eq. (9) may be approximated and expressed in the matrix 


form as follows: 


m=n 

W = [ E M{ e o} + E [ E(l n - < 10 > 

m= 1 

The continuum is next divided into small finite elements, and piece - 
wise continuous displacement fields are prescribed for each of such ele- 
ments in terms of their time -dependent nodal function values. Minimization 
of the total potential energy with respect to such parameters then yields the 
incremental equilibrium load -deflection equations of the entire structure. 
Such step-by-step incremental equations may finally be written in the global 
coordinate system: 

r n=n - 1 

K„-i]KH p 4 - E [ K n, ra -i]K4 -MW 

m=l 

n 

+ E (11) 

m= 1 
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with 



stiffness matrix derived from material matrix computed 

for the reduced time difference At-. = t ■ - t- 

b iJ b J 

external load vector at step n 
forces due to temperature changes 
body forces vector 


and in which the summation, as usual, signifies the memory of the material. 
The element stresses may then be obtained from Eq. (10), when element 
strains are derived from the usual relationship: 


— e 

U 


XU 


e 


— e 

u = au 
e = bu 

U , U being element nodal displacements in the global and local coordinate 
systems, respectively, \ the direction cosine matrix, and U and e are 
displacements and strains within the element. 

C. Special Incremental Procedure 

It is apparent from the nature of Eq. (11) that computation time may 
be excessive after a few time steps. This is because at each time step, 
recomputation of solution results is required for all preceding time steps, 
which are then added to obtain the final solution. However, in order to 
minimize such computation efforts, the program provides an option by which 
time steps may be so chosen that previous time intervals become a subset 
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of the following time intervals. Thus, the parameter £ may be expressed as 
the summation of incremental A£'s as follows (Ref. 2): 



i=l j=l 


where M defines the total number of time step groups, and N(i) is the number 
of steps in the ith group. The values of M and N can be suitably chosen by 
the user, and this scheme may be employed to solve the recursive Eqs. (10) 
and (11), provided material property and external loads are available for 
each £j. Figure 2 shows details of such a computation scheme for values of 
M = 3, N(l) = 2, N(2) = 3, and N(3) = 2. In such cases the time intervals 
tend to remain constant in the logarithmic scale; thus, it is then possible to 
cover a long time domain with relatively small computational effort. 

IV. INPUT PREPARATION FOR VISCEL 

The program VISCEL is assumed to be stored on a tape (say, 
number 12345), which contains the symbolic and relocatable program 
elements. Data deck corresponding to any problem must be preceded by 
a set of control cards which are first described below. The actual data deck 
preparation is explained next. 

A . VISCEL Control Cards for UNIVAC 1108/EXEC 8 Computer 

Depending on the size of the problem to be solved, the user may 
request for an appropriate core storage. Such values are assigned to an 
integer LDATA, to be calculated approximately from Ref. 2, Fig. 1, in 
which for most problems the major storage space would be required for 
elements of the upper symmetric half of the stiffness matrix. The program 
is compiled for LDATA = 20000 words decimal storage, and if more storage 
is requested, this is achieved by recompilation of two small programs 
COMBK and MAIN, the block data and the main driver programs, 
respectively. 

Furthermore, as explained in Ref. 2 (pp. 3-4), various Fastrand 
(drum) file storage units are utilized as additional stores during execution 
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of the program; their functions are summarized in Table A -7. Unless 
specified, the UNIVAC 1108 system automatically allocates 128 tracks to 
each of the units, which, however, may be inadequate for solution of large 
order problems. It is then necessary to increase the number of data tracks 
for such units by inserting relevant control cards in the run stream. 

Control cards corresponding to the two sets of values of LDATA are 
as follows: 

(1) Control cards with LDATA < 20000 

The following run stream may be used for problems which do 
not require more than 20000 words storage for the COMMON: 

@RUN,/TPC RUNID, ACCOUNT, PROJECT, TIME, PAGES 

@ MSG, READ TAPE 12345 

@ASG,T TAPE, T, 12345R 

@ FREE TPF$ 

@ASG, T TPF$, F///500 

@ COPY, G TAPE, TPF$ 

@ FREE TAPE 

|_@ASG,T_ UNIT NUMBER^/// 1000] 

@ XQT ABSEL 

YISCEL INPUT CARDS 
@ FIN 

(2) Card input with LDATA > 20000 (say, LDATA = 80000) 

When COMMON requirements are greater than 20000 (say, 
80000), the following typical run stream may be adopted: 

@RUN, /TPC RUNID, ACCOUNT, PROJECT, TIME, PAGES 

@ MSG, READ TAPE 12345 

@ASG,T TAPE, T, 12345R 

@ FREE TPF$ 

@ASG, T TPF$, F///500 
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@COPY, G TAPE, TPF$ 

@FREE TAPE 

[@ASCT T_ UNIT NUMBER7F2//7l000j 
@FOR, S COMBK, COMBK, COMBK 
- 2,2 

PARAMETER LDATA = 80000 
@FOR, S MAIN, MAIN, MAIN 
- 2,2 

PARAMETER LDATA = 80000 

@PACK 

@PREP 

@MAP, EN MAPEL, ABSEL 
@XQT ABSEL 
VISCEL INPUT CARDS 
@FIN 

Requests for additional storage tracks for the Fastrand units may be made 
by inserting control cards, shown above within the dotted boundaries. 

B. Input of Problem Data 

The physical arrangement of the data deck which follows the control 
cards (explained in previous section) is depicted in Fig. 3. This deck 
corresponds to values M = 2, N(l) = 2, N(2) = 3 in the time domain 
defined by Eq. (12). VISCEL input data may be as described below, with 
reference to Table 1 describing input items; the integers of the problem 
control card (Table 1, input item 2) is explained in Table A-8. 

Data Group 1 : Basic input for the elastic problem which also 

corresponds to the initial time solution of the 
viscoelastic problem 
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Data Group 2 ; Data for multiple solutions of the elastic problem 
or 

Data for viscoelastic incremental solution in the 
time domain 

The nature of the data in group 2, if any, is determined by the con- 
tents (ISUCA value) of the END card in the master (initial time) deck (input 
item 19, Table 1) and the subsequent additional input data decks for visco- 
elastic problems. Field specification for the END card is as follows: 

70X.I7, 3HEND (13) 

in which the 17 field corresponds to the integer ISUCA, which is to be set 
as follows : 


(1) 

ISUCA 

= 

0 

For linear elastic problems 

(2) 

ISUCA 

< 

0 

For multiple runs 

(3) 

ISUCA 

> 

0 

For linear viscoelastic problems 

(4) 

ISUCA 

= 

1 

For master and following deck 

(5) 

ISUCA 

> 

1 

For following additional decks in increasing 


sequence 


Thus, in the viscoelastic case, the first card following the data of 
the previous time step is the problem control card, equivalent to input 
item 2 of the initial time step, containing information on modifiable input 
items. The modified information is provided next, followed by the END 
card with the ISUCA value which determines the nature of the data, if any, 
in the succeeding step. A numerical example of a two-dimensional plane 
stress problem (Fig. 4) with irregular mesh labeling is chosen to elaborate 
on the preparation of the data; the complete input data are presented in 
Fig. 5 with M = 2, N(l) = 4, and N(2) = 2 values selected for the incre- 
mental time scheme of Eq. (12). Node relabeling may be requested by 
using appropriate option in input item 17 of Table 1. 

Relevant details on permanent and modifiable input items are provided 
in Table 2. Element data corresponding to input item 16 of Table 1 is 
described in Table A-9. Also, input items 13, 15, and 18 in the same table 
may be specifically described as follows: 
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Input item 13 (angle types - fixing local y and z axes) 


In connection with element type 4 (Table A-3), the input corre- 
sponding to column 1 6 of Table A-4 consists of a list of <j> angles 
in degree units. The <j> values are assigned quantities with abso- 
lute values less than 90 deg and are defined as the angle between 
the local y and global Y axes. Let the direction cosine vectors 

be denoted by (^ xX> ^ xY > ^ » ^ y X’ *yY’ * yZ ^ and ^zX’ * zY’ *zZ^ 
in which the local x axis is assumed to coincide with the nodal 

line 1-2 (Table A-3). Then the signs of -^ xX , ^ y x’ and ^zY are 
used to determine the sign of cj>; such procedure is summarized 
in Table A-10. 

Input item 15 (deflection boundary conditions) 

The deflection boundary condition relations may be written as 
(Ref. 1): 


“i,j = a 0 + a l U i',j’ + a 2V,j" + ••• (14) 

when coefficients a^, a ^ , a^, , and the input pairs (i,j), 

(i',j'), (i"»j")»...» are the relevant inputs as follows: 

i.j i»j a 0 

i.j i'.j' a 1 

i,j i",j" a 2 


in which the first two pairs along each row are the two degrees 
of freedom, under consideration and the related one, the third 
scaler relating such two deflection components. 
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(3) Input item 18 (concentrated load input) 

The inputs for the prescribed force boundary conditions for 
concentrated loads are as follows: 

i.j P 

i 1 » j ' P' 


where P. ., P* , , . , , are the prescribed concentrated nodal 
A # J 1 t j 

loads at nodes i, i * , . . . , corresponding to degrees of freedom 
j » j 1 »• • • . respectively. Apart from concentrated loads, the ele- 
ments may be subjected to any pressure as well as temperature 
loading as indicated in Table A-4. 

V. DESCRIPTION OF VISCEL OUTPUT 

Table A -11 provides a list of output items of the initial time step 
solution, whereas Table 3 summarizes such items for the entire visco- 
elastic problem with an input index value INP set to 1 for the elastic solution. 
The definition of stress components at mesh points is given in Table A -12. 

VI. ERROR MESSAGES AND DIAGNOSTICS 

The error messages shown in Table A-13 are usually related to the 
initial time step solution. Error message 10 in particular needs a detailed 
explanation, which appears either for geometrically unstable structures, or 
when the structure is not adequately supported. The last number appearing 
in the error message, if negative, indicates the mesh number to be checked 
carefully for existence of any unknown deformation. However, if the num- 
ber is positive, then it is first necessary to find from output item 10 of 
Table A -11 the pair of numbers with the second number identical to this 
error message number. The first number of the pair is called IBB, 
denoting the equation number in the reduced set of the stiffness matrix. 

Then column IBB of output item 10 of Table A-ll is searched for the row 
having the same IBB number found previously, such that the column IBO 
contains the number -1. The mesh number in that row happens to be the 
trouble spot, whereas the defective direction is the one appearing in the 
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table heading of the output item. In such case, the element descriptions, 
material matrices, and geometric continuity around the mesh point are to 
be checked to correct the situation. 

VII. CONCLUDING REMARKS 

This user's manual describes in detail the information necessary to 
utilize the computer program VISCEL for the solution of thermoviscoelastic 
problems associated with practical structures. Extensive applications of 
the problem are envisaged in the analysis of a wide variety of practical 
structures including solid propellant rocket motors, spacecraft components 
such as solar panels, etc. In order to make this document complete, some 
information, including most tables and figures in the Appendix, have been 
reproduced from Ref. 1. 
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Table 1. Input items (summary of options, contents, and formats) 


Input 

item 

No. 

Conditions 

determining 

options 

List of input statements that read the associated input item card(s) a 

Format 

(outside parentheses indicate the 
possibility of multiple cards) 

i 


(B±, i = 1, 14) The card may contain any alphanumeric message 

14A6 

2 


IN, IT, IDEG, ITYPE, IGEM, ISTR, IH, 18, IBN, IP, IPRS, IMAT, NTIC, 
ISDT, ISDY, ISDZ, IARE, IMMX, IMMY, IMMZ, IMFI, INX, INP, ISHUF, 
ICOR, IBUN, IMES, IPIR, ITAP, ITAS, Gl, G2, G3, ACEL 

214, 611, 314, 1012, 911, 3F5.4, E10.5 
(see Table A-8 for details) 

1 

ITYPE = 0 

(i, £;, G i , I = 1, IMAT) 

(3(12, 3E8.5)) 

ITYPE = 1 

(>, Djjj, Dj 2 ., D 14i< D 22 ., D 24 / D 4 4 .. D 55; D 56/ D 66/''“l/ 
a' ,i = 1, IMAT) 

(12, 9E8.5/I2, 2E8.5) 

ITYPE = 2 

(*/ Dll / £>12 , Dig , D i4 , ^15:' ^ 16 ;' D22 i '^23.' ^24 ^ 

til » i 1 

, D,„ , D,. , Do, , D,, , D. . , D,« , D., , i, D 5 , , D.. , 

26j' 33 ; ' 34^ 35/ 36/ 44/ 45/ 4e { ' 55/ 56/ 

D 60 , «i , a 2 ,a 3 .,i= 1,IMAT) 

(12, 9E8.5/I2, 9E8.5/I2, 6E8.5) 

4 

1 < IPRS < 99 

(f, p 4 , » = 1, IPRS) 

(8(12, E8.5)) 

IPRS = 0 

No input card 

5 

1 < NTIC < 99 

(i,h 4 ,i = 1, NTIC) 

(8(12, E8.5)) 

NTIC = 0 

No input card 

B 


(f, Alj, 1=1, ISDT) 

(8(12, E8.5)) 


No input card 

7 

9 

(i, (0l/0y)j, ; = 1. ISDY) 

(8(12, E8.5)) 

ISDY = 0 

No input card 

8 

1 < ISDZ < 99 

(/, (0F/0z)i. I = 1, ISDZ) 

(8(12, E8.5)) 

ISDZ = 0 

No input card 

9 

1 < IARE < 99 

( i,A if i = 1, IARE) 

(8(12, E8.5)) 

IARE = 0 

No input card 

10 

1 < IMMX < 99 

(/VC;,; = 1, IMMX) 

(8(12, E8.5)) 

IMMX = 0 

No input card 

11 

1 < IMMY < 99 

(/,/„.,/ = 1, IMMY) 
y i 

(8(12, E8.5)) 

IMMY = 0 

No input card 

12 

1 < IMMZ < 99 

(i,/„ i = 1 , IMMZ) 
i 

(8(12, E8.5)) 

IMMZ = 0 

No input card 

13 

1 < IMFI < 99 

(/, 0;,/' = 1, IMFI) 

(8(12, E8.5)) 

IMFI = 0 

No input card 

14 

ICOR = 0 
2 < IN < 9999 

X 

•< 

Nl 

II 

z 

(14, 3E12.4, 40X) or (40X, 14, 3E12.4) or 
(2(14, 3E12.4)) 

ICOR = 1 Input card(s) should be as required by the user's version of subroutine CORG (see Ref. 1) 
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Table 1 (contd) 


Input 

item 

No. 

Conditions 

determining 

options 

List of input statements that read the associated input item card($) n 

Format 

(outside parentheses indicate the 
possibility of multiple cards) 

15 

IBUN = 0 
1 < IBN < 9999 

(<k- h- ‘k- i'k- a ;.- k = ’' IBN ) 

(5(14,11, 14, 11, F6.0)) 

IBUN = 1 

Input card(s) should be as required by the user's version of subroutine BUNG (see Ref. 1) 

16 

IMES = 0 
1 < IT < 9999 

(MM m ,JlW,„,, J2W. m ,J3W m ,J4W. m ,J5W m m = 1, IT) 

(2014) (see Table A-9 for variables of 
the list) 

IMES = 1 

Input card(s) should be as required by the user's version of subroutine MESG (see Ref. 1) 

17 

ISHUF = Oor 1 

No input card 

ISHUF = 2 

(Nj, ; = 1 , in) 

(2014) 


(Nj, /MAX;, i = 1, IN) 

(2014) 

18 

1 < IP < 9999 

(/(,/,, P,,l= 1, IP) 

(5(14, 11, El 1.4)) 

ir = 0 

No input card 

19 


No list (the card is punched END in the last three columns) 

70X, 17, 3HEND (ISUCA value) 

20 


No input for standard VISCEL; otherwise input of certain user's subroutines (see Ref. 1) 

VISCEL PROBLEM CONTROL CARD, PROVIDES MODIFIABLE INFORMATION 


MODIFIED INFORMATIONS 

I 


I 


19 


END card 

70X, 17, 3HEND (ISUCA value) 


PROCESS TO BE REPEATED 




FOR SUBSEQUENT TIME 

STEPS 


a Nomenclature 




P 

pressure 

f s 

moment of inertia about local z axis 


h 

thickness 

<t> 

angle determining the orientation of principal axes of cross 

section 

A t 

temperature increase 


in overall coordinate system 


d t/dy 

temperature gradient in local y-axis direction 

X, Y, Z 

overall coordinates of mesh points 


6 f/d z 

temperature gradient in local z-axis direction 

x, y, * 

local coordinates 


A 

cross-sectional area 

I'*,/*), I'*',/*')' a k 

index pairs and the constant of the kth dbc input unit (see 
111 -D) 

Section 

C 

torsional constant 

'riiSi 

index pair and constant of the / th concentrated load input unit (see 

l t/ 

moment of inertia about local y axis 


Section IV-B) 



b The symbols shown in Input Item 3 are defined in Figs. A-2c, 2d, and 2e. 
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Table 2. Permanent and modifiable input items 


Input item 
number 

Description of input item 

Qualifications 

Existence 
in the 
master 
deck 

Existence 
in the 
successive 
decks 

1 

Title card 

For master deck only 

Jje 

- 

2 

Control card 

For master and successive decks 


9j< 


Modifiable information 


3 

Material types 

For master and successive decks 

* 

o 

4 

Pressure types 

For master and successive decks 

O 

o 

5 

Thickness types 

For master and successive decks 

O 

o 

6 

Temperature increase types 

For master and successive decks 

o 

o 

7 

Temperature gradient types 
-local y-axis direction 

For master and successive decks 

o 

o 

8 

Temperature gradient types 
-local z-axis direction 

For master and successive decks 

o 

o 

9 

Cross-sectional area types 

For master and successive decks 

o 

o 

10 

Torsional constant types 

For master and successive decks 

o 

o 

11 

y-moment-of-inertia types 

For master and successive decks 

o 

o 

12 

z-moment-of-inertia types 

For master and successive decks 

o 

o 

13 

Angle types-fixing local y and 
z axes 

For master and successive decks 

o 

o 

Permanent information 



■■ 

Mesh point coordinates 

For master deck only 

o 

- 

PS ^ 

Deflection boundary conditions 

For master deck only 

o 

- 

in 

Element descriptions 

For master deck only 

o 

- 

17 

Relabelling information 

For master deck only 

o 



Modifiable information 


18 

Concentrated loads 

For master and successive decks 

o 

o 

19 

End card 

For master and successive decks 

* 


* = the card(s) must exist. 




o = the card(s) exist optionally depending on 
modifiable information. 

the control constant in the control card and relate 

to the 

- = the card(s) must not exist. 
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Table 3. Summary of output items 


Output 

item 

number 

Output for linear elastic solution or 
initial time solution of linear 
viscoelastic problems 

Output for linear 
viscoelastic solution 

1 

Linear elastic problem or linear visco- 
elastic problem 

(1) Title of the problem 

(2) Table for control constants 

Number of equal time step 
group 

Number of time steps in the 
group. . . . 

2 

Modifiable information (material prop- 
erties, pressure types, etc.) 

Modifiable information (mate- 
rial properties, pressure 
types , etc. ) 

3 

Nodal coordinates 


4 

Mesh topology; element property types 


5 

Relabelling message 


6 

Topology of the reduced stiffness matrix 


7 

Stiffness matrix requires. . . . storage 
locations 


8 

Total common length is (decimal). . . . 
storage locations 


9 

Count of main diagonal elements of row 
listed stiffness matrix 


10 

Force and displacement boundary condi- 
tions indirections 1 (2, 3, 4, 5, 6) 


13 

Input link took. . . . seconds 

Input link took. . . . seconds 

17 

Generation link took. . . . seconds 

Generation link took. . . . seconds 

19 

Nodal deflections 

Accumulative nodal deflections 

20 

Forces acting at the nodes 

Forces acting at the nodes 

21 

Deflection link took. . . . seconds 

Deflection link took. . . . seconds 

22 

Stresses at the nodes 

Stresses at the nodes 

25 

Stress link took. . . . seconds 

Stress link took. . . . seconds 
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Fig. 3. 


Physical arrangement of data deck for the 


VISCEL program 
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Fig. 5 (contd) 
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Table A-l. Deflection degrees of freedom at a point 
for different cases of structures 


Column number 

F 

1 2 

ID 

ID 

D 

6 

F 

Case number 

\ , E = 

\ 2 -o o 

\ k 9 £ 

Case \ 

description \ 

Displacement along X 

• 

Displacement along Y 

Displacement along Z 

Rotation about X 

Rotation about Y 

Rotation about Z 

Number of degrees of freedom 

i 

Planar truss 

11 

HI 



■ 


2 

2 

Space truss 

89 

■ 

§j 


D 


3 

3 

Planar frame 

H 

■ 

II 


D 


3 

4 

Space frame 

■ 

■ 



■ 


6 

5 

Gridwork frame 







3 

6 

Plane stress 

■ 

■ 

D 


D 


2 

7 

Plane strain 


■ 





2 

8 

Plate bending 



■ 


■ 


3 

9 

General solid 

111 

■ 



D 


3 

D 

General shell; bend., memb. 

111 


1 


U 


6 

n 

General shell, membrane 

fjj 

■ 

jjj 


u 


3 

,2 

Solid of revolution 

jj| 

U 

D 




2 


Shell of revolution, membrane 

jj§ 


■1 


m 


2 


Shell of rev.; bend., memb. 

m 

■ 

□ 


J 


3 


°X r Y, Z refer to the axes of the overall coordinate system. 
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Case number 


Table A-2. Types of structures that VISCEL can handle 
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Shell of revolution, membrane 
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OJ 
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Table A-3. 


DI 

2 

3 

4 

5 

6 

7 


Element type number 

Element geometry 

Number of nodes (vertices), IMS 

Degrees of freedom per node, IDEG 

Number of words for element 
description, 18 

Case No. of structure (Table 111-2) 
for which this element may be used 

Nodal line on the first material 
axis direction 


1 

Line segment 

2 

2 

5 

1 

1-2 


1 

Line segment 

2 

3 

5 

2 

1-2 


2 

Line segment 

2 

3 

6 

3 

1-2 


3 

Line segment 

2 

3 

5 

5 

1-2 


4 

Line segment 

2 

6 

8 

4 

1-2 


5 

Triangle 

3 

2 

6 

6,7 

• 


6 

Quadrilateral 

4 

2 

a 

mm 

• 


7 

Triangle 

3 

3 

6 

8 

• 


8 

Quadrilateral 

4 

3 

7 

8 

• 


9 

Tetrahedron 

4 

3 

6 

9 

• 


10 

Hexahedron 

8 

3 

10 

9 

• 


11 

Triangle 

3 

6 

6 

10 

1-2 


12 

Quadrilateral 

4 

6 

7 

10 

1-2 


13 

Triangle 

3 

3 

6 

11 

1-2 


14 

Quadrilateral 

4 

3 

7 

11 

1-2 


15 

Triangular torus 

3 

2 

6 

12 

• 


16 

Quadrilateral torus 

4 

2 

7 

12 

• 


17 

Conical segment 

2 

2 

5 

13 

1-2 


18 

Conical segment 

2 

3 

5 

14 


1 


Element properties 



Legend 


structure is in the overall (X-Y) plane 

the mesh is in the overall (X-Y) plane and overall 
Y axis is the axis of revolution 

first material axis is the overall X axis 

normal to nodal line 1-2 and the overall Z axis, 
and away from element 

normal to surface shown in column 8 and in direc- 
tion of local normal 

local z-axis direction 

Perpendicular to the element in the plane estab- 
lished by the element and the overall X axis. The 
direction is such that the angle between the per- 
pendicular and the X axis is less than 90 deg 


1-2-3 


1 -2-3-4 


1-2-3 


1-2 -3 -4 


1-2 


1-2 


1-2 (*) 


1-2 (*) 


Any 


Any 


Any 


Any 


© 


© 


© 


© 


A 


A 


A 


A 


▲ 


A 


f in the direction of overall Z axis 

■ parallel to overall axes (for element type 1, for 
stresses, local system as in □) 

A x axis: nodal line 1-2; z axis: normal to middle 
surface, which sees labels counterclockwise 

A x axis: nodal line 1-2; y axis parallel and opposite 
to Z axis 

D x axis: nodal line 1-2; y axis: one of the principal 
axes of the cross sections 


































































































































































Table A -4. Necessary and optional information for element definition 
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Temperature gradient along y type number 

















Table A-5. Types of elements available for 
different cases of structures 


f\ 

Column number 

III 

ID 

D 

D 

ID 

D 

ID 

IB 

Case number 

\ a. 

N. X 

N. C 

\ © 

\ E 

Case N. m 

description 

Line segment 

Triangle 

1 


Tetrahedron 

Hexahedron 

D 

0 

JO 

3 

0> 

c 

JO 

h- 

Quadrilateral torus 

i 

Planar truss 

1 








2 

Space truss 

i 






B 


3 

Planar frame 

g 

D 



D 


B 


4 

Space frame 

4 








5 

Gridwork frame 

3 








6 

Plane stress 


a 

g 






7 

Plane strain 


D 

m 






8 

Plate bending 


§3 

m 

m 

D 




9 

General solid 


i 

■ 

u 

B 




10 

General shell; bending, membrane 


■ 

B 

i 

B 

B 



11 

General shell, membrane 


13 




B 



12 

Solid of revolution 


i 

i 




a 

B 

13 

Shell of revolution, membrane 




n 



m 

B 

14 

Shell of revolution; bending, membrane 




B 

1 


□ 

□ 
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Table A-6. Convention for ordering the 
vertices of elements 


Element 

type 

number 

First 

vertex 

Other vertices 

i 

Any 

The remaining 

2 

Any 

The remaining 

3 

Any 

The remaining 

4 

Any 

The remaining 

5 

Any 

Counterclockwise sequence about overall Z axis 

6 

Any 

Counterclockwise sequence about overall Z axis 

7 

Any 

Counterclockwise sequence about overall Z axis 

8 

Any 

Counterclockwise sequence about overall Z axis 

9 

Any 

Counterclockwise sequence for the first three vertices 



about the normal of their plane, heading towards 
the fourth vertex 

10 

Any 

* 

11 

Any 

Clockwise sequence about local normal** 

12 

Any 

Counterclockwise sequence about local normal** 

13 

Any 

Counterclockwise sequence about local normal** 

14 

Any 

Counterclockwise sequence about local normal** 

15 

Any 

Counterclockwise sequence about overall Z axis 

16 

Any 

Counterclockwise sequence about overall Z axis 

17 

* * * 

The remaining 

18 

*** 

The remaining 

^Counterclockwise sequence for the first four vertices on the same face 

about the normal heading towards the other four vertices. The fifth 
vertex lies diagonally across the first vertex. The last four vertices also 
establish a counterclockwise sequence about the normal of their face, 
heading towards the first four vertices. 

! ** Local 

normals 

lead always to the same side of the space divided by I 

the middle surface. 

***The one with smaller meridional arc length (the meridional curve should 

have 

a direction). 
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Table A-7. The functions of the FORTRAN units as used in VISCEL 


FORTRAN 

unit 

number 

Function of the unit 

1 

System 

2 

Chain 

3 

Scratch for topological information generated in Link 3 

4 

Storage for deflections 

5 

Input 

6 

Output 

7 

Punch 

8 

Overlays for FORTRAN IV 

9 

Storage for material elastic constants 

10 

Storage for material expansion coefficients 

11 

Storage for temperature changes 

12 

Storage for elemental stiffness matrices 

13 

Storage for overall stiffness matrices 

14 

Scratch for elemental stiffness matrices 


Scratch for incremental deflections for ISTEP = 1 

15 

Storage for stiffness matrix decomposed by Choleski scheme 
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Table A-8. Summary of the problem control card (input item 2) of input data 


Name of 
field 

Card columns 
of field 

Format 

Range 

Description 

IN 

1-4 

14 

2-9999 

Total number of mesh points 

IT 

5-8 

14 

1-9999 

Total number of elements 

IDEG 

9 

11 

2-6 

Number of degrees at a mesh point (see Table A-l, column 7) 

ITYPE 

10 

11 

0-2 

Material indicator: 0 — isotropic, 1 — orthotropic, 2 — general \(see Fig. A-2) 

IGEM 

11 

11 

0-1 

Geometry indicator: IGEM = 0 all Z coordinates are zero" (see Table A-3, 
IGEM = 1 not all Z coordinates are zero" column 10) 

ISTR 

12 

11 

0-1 

Plane strain case indicator: ISTR = 1 plane strain 

ISTR = 0 not plane strain 

IH 

13 

11 

2-8 

Maximum number of vertices in elements used (see Table A-3, column 3) 

18 

14 

11 

5-1 0 b 

Maximum number of words for element description ((see Table A-3, column 5) 

IBN 

15-18 

14 

1-9999 

Total number of dbc input units (see Section III) 

IP 

19-22 

14 

0-9999 

Total number of concentrated load input units (see Section III) 

IPRS 

23-26 

14 

0-99 

Total number of different pressures 

IMAT 

27-28 

12 

1-99 

Total number of different materials 

NTIC 

29-30 

12 

0-99 

Total number of different thicknesses 

ISDT 

31-32 

12 

0-99 

Total number of different temperature increases 

ISDY 

33-34 

12 

0-99 

Total number of different temperature gradients ©f/0 y)' 

ISDZ 

35-36 

12 

0-99 

Total number of different temperature gradients (0f/0z)' 

IARE 

37-38 

12 

0-99 

Total number of different cross-sectional areas 

IMMX 

39-40 

12 

0-99 

Total number of different torsional constants 

IMMY 

41-42 

12 

0-99 

Total number of different moments of inertia (about y axis) c 

IMMZ 

43-44 

12 

0-99 

Total number of different moments of inertia (about z axis)' 

IMFI 

45-46 

12 

0-99 

Total number of angles fixing local y and z axes c 

INX 

47 

11 

1-4 

Number of link after which return-to-beginning-for-next-job is done 

INP 

48 

11 

0-2 

Printout indicator: 0— minimum; 1 —intermediate; 2— detailed output 
(see Table A-l 1) 

ISHUF 

49 

11 

0-3 

Relabelling indicator: 0— no relabelling; 1— iterate to relabel without reading 

cards; 2— read cards and iterate to relabel; 3— relabel 
as shown on cards (see Ref. 1) 

ICOR 

50 

11 

0-1 

Indicator for coordinate generation: 0— read coordinates from cards; 

1— generate coordinates via subroutine CORG (user's version) (see Ref. 1) 

IBUN 

51 

11 

0-1 

Indicator for displacement boundary conditions: 0— read from cards; 
1— generate with user's version of subroutine BUNG (see Ref. 1) 

IMES 

52 

11 

0-1 

Indicator for element descriptions: 0— read from cards; 1— generate with user's 
version of subroutine MESG 

1 PI R 

53 

11 

0-2 

Local coordinate selection indicator for shells: 0— assume local x as 1-2 line 
of lowest numbered element; 1— assume as principal; 2— read by subroutine 
AGEL 
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Table A-8 (contd) 


Name of 
field 

Card columns 
of field 

Range 

Format 

Description 

ITAP 

54 

ii 

0-9 

Chain tape number for program (if zero, program assumes 2) 

HAS 

55 

ii 

0-9 

Chain tape number for intermediate storage 

G1 

56-60 

F5.4 

(-1.H+1.) 

Cosine of the angle of acceleration vector with X axis" 

G2 

61-65 

F5.4 

(-1-H+1.) 

Cosine of the angle of acceleration vector with Y axis" 

G3 

66-70 

F5.4 

(— 1-M+l.) 

Cosine of the angle of acceleration vector with Z axis" 

ACEl/ 1 

71-80 

El 0.3 

Any 

Magnitude of acceleration vector times unit mass (unit weight) 

S, X, Y, Z refer to overall coordinate system. 



b When 18 = 10, zero should be punched in column 14. 


c x, y, z refer to the local coordinate system of the element. 


d ln element type 3, ACEL means weight per unit length. 
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Table A-9. Description of element data for different element types 


Element 

Meanings of variables appearing in the list of Input Item 16 of Table IV-1“ 


type No. 

MM 

J1W 

J2W 

J3W 

J4W 

J5W 

J6W 

J7W 

J8W 

J9W 

J10W 

1 

M 

1004- IMET 

100(JARE) + ITEM 

JPRS 

N, 

N 2 






2 

M 

200 + (MET 

100(JARE) + ITEM 

JPRS 

100(JMMZ) + JSDY 

Ni 

n 2 





3 

M 

300+ IMET 

100(JPRS) + JSDZ 

100(JMMX) + JMMY 

Nt 

N, 






4 

M 

400 + IMET 

100{JARE) + ITEM 

1 00(J MMX) + JMMY 

1 00 ( J MMZ) + J SDY 

100(JSDZ)+JMFI 

JPRS 

N, 

n 2 



5 

M 

500 + IMET 

100(ITIC)+ ITEM 

JPRS 

Ni 

n 2 

N s 





6 

M 

600+ IMET 

100(ITIC) + ITEM 

JPRS 

N, 

n 2 

N 2 

n 4 




7 

M 

700+ IMET 

T 00(ITIC) + ITEM 

100(JSDZ) + JPRS 

Ni 

N 2 

n 3 





8 

M 

800 + IMET 

1 00(ITIC) + ITEM 

100(JSDZ) + JPRS 

Ni 

n 2 

n 2 

n 4 




9 

M 

900+ IMET 

100(JPRS) + ITEM 

N, 

N 2 

N, 

N, 





10 

M 

1000+ IMET 

100(JPRS)+ ITEM 

N, 

n 2 

N a 

n 4 

Ns 

N« 

N, 

Nh 

11 

M 

1 100+ IMET 


100(JSDZ) + JPRS 

Ni 

n 2 

n 3 





12 

M 

1200 + IMET 


100(JSDZ) + JPRS 

N, 

N, 

n 2 

n 4 




13 

M 

1300+ IMET 

1 00(IT1C) + ITEM 

JPRS 

N, 

N 2 

n 3 





14 

M 

1400 + IMET 

100(ITIC) + ITEM 

JPRS 

N, 

n 2 

n 3 

N ( 




15 

M 

1500+ IMET 

1 00(ITIC) + ITEM 

JPRS 

N, 

n 2 

n 3 





16 

M 

1600+ IMET 

1 00(ITIC) + ITEM 

JPRS 

N, 

N, 

n 3 

N, 




17 

M 

1700+ IMET 

100(ITIC)+ ITEM 

JPRS 

N, 

N, 






18 

M 

1800+ IMET 

1 OO(ITIC) + ITEM 

100(JSDZ) + JPRS 

N, 

n 2 





— 


“Definitions: 

IMET material type number 
ITIC thickness type number 
JPRS pressure type number 
JARE cross-sectional area type number 
ITEM temperature increase type number 


JSDY temperature gradient along y type number 
JSDZ temperature gradient along z type number 
JMMX torsion constant type number 
JMMY moment of inertia (about y) type number 
JMMZ moment of inertia (about z) type number 


JMFI angle fixing local y, z axes type number 
Nj mesh-point label of the first vertex 
N 2 mesh-point label of the second vertex, . . . 
N g mesh-point label of the eighth vertex 
m element label 


M = — ( j 000 ) * which is to be interpreted in FORTRAN integer arithmetic sense 

Note: On the cards of Input Item 16 (first option), zero field(s) (of four columns) after the last nonzero field of 
element description is ignored if provided. 


-J 
































































































































Table A -10. Table for determining the direction of local y axis 

and the sign of angle $ 
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Output 
item No. 


Table A-ll. List of output items 


Output 
item No. 

Output item 

Minimum 
INP = 0 

Intermediate 
INP = 1 

Detailed 
INP = 2 

Link 1 (input link) 

i 

Table for title and important constants 

■ 

■ 

■ 

2 

Tables for material, loading, and geometry 


11 

1 

3 

Table for coordinates of mesh points 


■ 

§Jj 

4 

Table for element properties 


■ 

■ 

5 

Message and/or tables and punched cards 
for relabelling 

|§§ 


@ 

6 

Table for reduced stiffness matrix 


1 

■ 

■ 

Message about necessary storage for re- 
duced stiffness matrix 

11 

I 

■ 

8 

Message about total common length 


■ 


9 

Table for diagonal elements of reduced 
stiffness matrix 

■ 

n 

m 

10 

Tables for force and deflection boundary 
conditions 

■ 

u 

■ 

11 

Table for common, integer (IA block) 



■ 

12 

Table for common, floating (AA block) 

■ 

■ 

■ 

13 

Message for the execution time of Link 1 


■ 

■ 


Output produced (blank means no output produced). 

@ Output related with relabeling 

* May be produced selectively by subroutine CAS2 (see Ref. 1). 
t May be produced selectively by subroutine CAS4 (see Ref. 1) 


Output 
item No. 

Output item 

Minimum 
INP = 0 

Intermediate 
INP = 1 

Detailed 
INP = 2 


Link 2 (generation link) 

14 

Tables for element stiffness matrices 

* 


- 

15 

Table for upper half of reduced 
stiffness matrix 

* 

■ 

■ 

16 

Table for reduced load vector 

* 

D 

i 

17 

Message for the execution time of Link 2 


^8 

: 

■ 

Link 3 (deflection link) 

18 

Table for reduced solution vector 




19 

Table for deflections 

: ' 

MB 

||i 

20 

Table for forces acting on the nodes 


■ 


21 

Message for the execution time of Link 3 


■ 

I./ 

■ 

Link 4 (stress link) 

22 

Table for stresses at the nodes 

m 

■ 

; 

23 

Details of the best-fit stress computation 

D 

a 


24 

Table for stresses of one-dimensional 
elements 



u 

25 

Message for the execution time of Link 4 
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Intermediate 

























































Table A-12. Meanings of the components of stresses at mesh points of 
two- and three-dimensional continua 


Class 

No. 

Element 

type 

No. 

Structure 

type 

First 

component 

Second 

component 

Third 

component 

Fourth 

component 

Fifth 

component 

Sixth 

component 

i 

5,6 

2-D elasticity 


02 

r 12 * 

o* 

0 

0 

2 

7,8 

Plate, bending 

0 

0 

0 

Mi 

m 2 

Mu 

3 

15, 16 

Solid of revolution 

Ol 

02 

o 3 

T 12 

0 

0 

4 

9, 10 

General solid 

Oi 

02 

0 3 

T 12 

Tl3 

T23 

5 

17 

Shell of revolution, membrane 

N, 

n 2 

0 

0 

0 

0 

6 

18 

Shell of revolution, membrane and 
bending 

N, 

n 2 

0 

M, 

M : 

0 

7 

13, 14 

Shell, membrane 

N, 

n 2 

Nj2 

0 

0 

0 

8 

11, 12 

Shell, membrane and bending 

N, 

n 2 

Nj 2 

M, 

m 2 

Mu 

“Nomenclature: 

°i> °s normal stresses N y N 2 membrane normal forces M y M 2 bending stress couples 

T i 2 ' T i3' T 23 s ^ ear stresses N 12 membrane shear M J2 twist stress couple 

*CTj for 
t TlI for 

plane-strain cas 
plane-strain ca 

e 

le 


1 

% 

r 

rJ 

T 2y / 

71 

23 7 

* ► 2 . 

v / y 


2 

1 

E 

OUPLES 

/ 

1 

The axis labe 

STRE 

s 1, 2, a 

K X N n X /If, 

SS 1 MEMBRANE FORCES 1 STRESS C 

nd 3 stand for KSI, ETA and ZTA {direction cosines for local axes), respectively 
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Table A-13. List of error messages 


No. 

Error message 

1. 

INPUT ERROR 

2. 

THE FOLLOWING DISPLACEMENT BOUNDARY CONDITION(S) 
CAUSE(S) MORE THAN ONE MULTIPLE CONNECTION FOR 
THE UNKNOWNS. THEY ARE IGNORED 

3. 

i: IN ELEMENT .... ERROR IN MESH TOPOLOGY INFORMA- 
TION. NO CORRECTION IS MADE, if: IN ELEMENT ....... 

PROPERTY TYPE NUMBER(S) IS OUTSIDE THE RANGE. THE 
TYPE NUMBER(S) IS ASSUMED LARGEST POSSIBLE 

4. 

ELEMENT ... IS UNACCEPTABLE. DISREGARDED 

5. 

WARNING. LESS THAN 12750 DECIMAL LOCATIONS ARE 
AVAILABLE FOR THE NEXT LINK PROGRAMS. THOUGH IT 
MAY BE SUICIDAL, EXECUTION CONTINUES 

6. 

THE POINT . . . DOES NOT APPEAR IN THE MESH 
TOPOLOGY 

7. 

DUMMY AREA OVERLAPS COMMON AREA BY . . . DECIMAL 
LOCATIONS. RECOMPILE BY CHANGING THE EQUIVA- 
LENCES OF DUMMY AND BB IN LINKS 1 AND 3, RE- 
SPECTIVELY 

8. 

ELEMENT IS UNACCEPTABLE. DISREGARDED . . . 

9. 

THE VOLUME OF ELEMENT IS TOO SMALL . . . 

DISREGARDED 

10. 

STIFFNESS MATRIX IS NOT POSITIVE DEFINITE . . . 

11. 

NO SCRATCH TAPE IS GIVEN OR ERROR IN SCRATCH TAPE 

12. 

MORE THAN 12 NON-ONE-DIMENSIONAL ELEMENTS AT 
NODE . . . 

13. 

NODAL STRESS COMPUTATION IS DELETED DUE TO 
PRECEDING 

14. 

NO SCRATCH TAPE. STRESS LINK IS NOT EXECUTED 

15. 

ERROR IN READING ELEMENT SETS FROM TAPE ITAS. STRESS 
LINK EXECUTION IS DELETED 

16. 

NOT ENOUGH INDEPENDENT INFORMATION AVAILABLE 

17. 

ERROR IN MESH TOPOLOGY. NODE ASSUMED INTERNAL 

18. 

MORE THAN 4 MATERIALS, FIRST 4 CONSIDERED 

19. 

MORE THAN 4 CLASSES, FIRST 4 CONSIDERED 

20. 

MORE THAN 19 ELEMENTS, FIRST 19 CONSIDERED 

21. 

NOT ENOUGH INFORMATION FOR BEST-FIT QUADRATIC 
BEST-FIT PLANE IS USED 

22. 

NOT ENOUGH INFORMATION FOR MIDDLE SURFACE NOR- 
MAL. APPROXIMATE XII AND ZTA VALUES ARE USED 

23. 

SCRATCH AREA FF OVERLAPS WITH RESIDUAL AREA. PUSH 
FF FURTHER DOWN BY RECOMPILING LINK 4. 
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ISOTROPIC CASE 


[D,] = 


aE bE bE 0 0 0 

aE bE 0 0 0 

aE 0 0 0 

GOO 
SYM G 0 

G 


G (4G — E) G (2G) 

, b = 


E (3G-E) 
[a] = [or, a, a] 
INPUT: E, G, a 


E (3G-E) 




ORTHOTROPIC 

CASE 





£>;, d;, * 

DU 

0 

o' 




DU * 

DU 

0 

0 


r d 

i = 

t 

0 

0 

0 


L 

J 


d;, 

0 

0 




SYM 


Dt 5 

D£g 




- 



DU 


* 

IN PLANE-STRAIN CASE 

Dn, OTHERWISE 0 

t 

IN PLANE-STRAIN CASE 

NOT NECESSARY, 


OTHERWISE 0 





[a] = 

[arl, ai, 0] 





INPUT 

DU, D'i 2, Dt 4, DU, 

DU, 

DU, 

Dm, D(jn 



at, ai 





GENERAL CASE 


Dll Di 2 D 13 Dl 4 Dis 

Die 



D22 D23 D24 D 25 

D26 



D 33 D34 D33 

Dja 


[ D 0 = 

D44 Dis 

Dio 



SYM D55 

D50 




Dee 






[a] = 

[an, ai, ar 3 ] 


INPUT: 

D 12 , Di 3, D 14, Dl5/ Did, D- 22 , 


? 23 / D24/ D25/ D2G# D33, D34, D35, 


^ 36 / D44/ Dis, D4G, D 5 3, D 5C , £ 

6G 

ott, ai, «3 



(c) 


(d) 


(e) 


Fig. A-2. Description of the material 
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